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Abstract 

The purpose of these notes is to give some examples illustrating how naive numerical approx- 
imations to PDE's may not work at all as expected. In addition, the following two important 
notions are introduced: (I) von Neumann stability analysis — helps identify when (and 
if) numerical schemes behave properly. (II) Artificial viscosity — a tool in stabilizing nu- 
merical schemes. These notes should be read in conjunction with the use of the MatLab 
scripts (in the Athena 18311-Toolkit at MIT) whose names end with the acronym GBNS (for 
Good-Bad-Numerical-Schemes) . 
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1 Naive Scheme for the Wave Equation. 

We will illustrate the points we want to make with the wave equation (in one space dimension) 

92U ^ 0. (1.1) 



dt 2 dx 2 

Since this equation is second order in time, it needs two initial conditions. For example: 

du 

u(x, 0) = uo(x) and — (x,0)=vo(x). (1.2) 
We will assume here that both u and v are periodic, with some period | T > 0. | Then the solution 



of (1.1) is periodic in x with the same period: u(x + T,t) = u(x, t). 
Remark 1.1 We note that, in fact, we can write the solution of this problem explicitly 

U= \ ~~ ^ + u °( x + ^ + / v o(s)ds S j . 

However, this is not the point here (see below). 

Operate now as if (1.1) were complicated enough that we needed to solve the equation numerically. 
For this purpose introduce a numerical grid {x n ,tj} — where n and j are integers, as follows 

x n = xo + nAx and tj = jAt . (1.3) 

Here Ax and At are some "small" positive constants and x$ is arbitrary. Next replace the function 
u = u(x,t) of the continuum variables x and t by a discrete double sequence {u J n }, where 

u 3 n = u{x n ,t j ) . (1.4) 
Finally, introduce the new variable v = — to re- write equation (1.1) as a first order in time system 

du dv d 2 u 

m =v and m = d^- (L5) 

In view of (1.4) it is now clear that u J n (and the similarly defined v 3 n ) should satisfy 

n M n =vi + 0(At) and = ^ +1 ( ^ )2 + ^ + 0(At, (Ax) 2 ) , (1.6) 

which can be checked by expanding u 3 n+1 , ... in Taylor series centered at (x n , tj) — using (1.4) 
— and substituting the expansions in (1.6). This suggests the following numerical scheme, allowing 
simple calculation of the solution at time t = tj + i (once it is known at time t = tj) 

At_ 
(AJ) 

where the errors should be of size 0(At, (Ax) 2 ), that is: small. 

Upon implementation one quickly discovers that this algorithm is disastrously bad. The MatLab 
scripts: InitGBNS, lectureGBNS, demoGBNS, movieGBNS and the help file readmeGBNS in the Athena 
18311-Toolkit all deal with this scheme and another one to be introduced later in these notes. In 
particular, lectureGBNS goes through and explains a series of calculations showing the details of 
how the scheme fails. We illustrate here the problem with a couple of examples. 



< +1 = < + At < and < +1 = 4 + -— - 2< + , (1.7) 
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Example 1.1 Consider the following initial data (with period T = 2) for equation (1.5): 



u(x, 0) = Uq(x) = - (1 + cos(7T2;)) and v(x, 0) = vq(x) = 

At 



The exact solution: 



U = - (2 + C0S(7T (x — t)) + COs(7r(x + t))) = - (1 + C0S(7T X) C0S(7T t)) 



(1.8) 



see 



remark 1.1 — is clearly also periodic in time of period 2 (a standing wave). For the numerical 
solution we take 



Ax = 2 At = 2/N (for some "large" N) and 



x 



1 1 in (1.3). Then we im- 



are equivalent) and solve the equations over one time period: 



plement (1.1) for 1 < n < N (the periodicity of the solution means that the indexes n + N and n 

< t < 2. | 




Figure 1.1: Solution of (1.5) with initial data (1.8) using (1.7) with 40 points in the 
space grid. To avoid an over-dense graph not all the points in the numerical grid are 
plotted. However, enough points to show all the relevant details are kept. 



Figure 1.1 shows the result of this calculation using N = 40. Note that the periodicity in time fails 
to hold. In fact, after one time period the numerical method appears to have amplified the initial 
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data by about 30%/ However, maybe this is not so bad (or is it?); after all the value of N being 
used is not that large and the numerical solution looks otherwise quite reasonable. 

Let us now check what happens as we increase the resolution (larger N). Any reasonable numerical 
scheme ought to give a better approximation when we do this. Figure 1.2 shows the result of in- 
creasing N to N = 57 (a rather small increase). The new approximation is not only not better; it 
is a disaster. By time t ~ 2, 0(1) grid scale (i.e. wavelength = 2 Ax) oscillations appear in the 
numerical solution, making it useless. As we will soon see, the scheme is amplifying the errors; the 
30% amplification of the initial cosine wave seen when using N = 40 was just a forewarning of what 
happens for larger N. As N is made even larger, the oscillations generated become huge (in fact, 




Timet — dt=1/N. Space x — dx=2/N. 

Figure 1.2: Solution of (1.5) with initial data (1.8) using (1.7) with 57 points in the 
space grid. To avoid an over-dense graph not all the points in the numerical grid are 
plotted. However, enough points to show all the relevant details are kept. 



their size increases exponentially with N, as we will soon show). This is illustrated by figure 1.3, 
which corresponds to N = 80. Here (instead of a 3D graph) we plot the numerical solution at time 
t = 2. Grid scale (wavelength = 2Ax) oscillations is all that can be seen in this graph — notice the 
(very large) vertical scale on this figure! 
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Finally, we point out that if (instead of increasing N ) we compute for longer times, the same effect 
of large amplitude grid scale oscillations arising (which grow exponentially in time) is observed. 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



Space x — dx=2/N. Solution for time t = 2 

Figure 1.3: Solution of (1.5) with initial data (1.8) using (1.7) with 80 points 
in the space grid. Notice the large amplitude grid scale oscillations generated by 
the scheme. There is nothing but numerical noise in this picture! 



Example 1.2 In a second example we take the following Gaussian initial data for equation (1-5) 

u(x, 0) = u (x) = exp(— a ln(10) x 2 ) and v(x, 0) = vq{x) = , (1.9) 

for — 1 < x < 1, where a > is a constant. We extend this to periodic initial data (of period T = 2) 
by repeating the above profiles over each interval (2n — 1) < x < (2n + 1), with n integer. These 
initial values are not smooth — as were the ones in the prior example. There is a small corner in 
uq{x), whenever x is an odd integer (in particular for x = ±1 ). This is because at these points there 
is a cut-off from a Gaussian centered at x — 1 to one centered at x + 1. Notice that the size of the 
miss-match in the derivatives of u goes down very rapidly as a increases. 
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For the numerical solution we take x$ = —1, Ax = 0.02 and At = 0.01 in (1.3) — this corresponds 
to N = 100 in the notation of example 1.1 — and use (1.7) to solve the equations for < t < 0.5. 
This is very similar to what we did in the prior example, except that here we vary the initial 
conditions (by changing the parameter a) instead of changing the resolution with variations in N . 

In the first calculation, we take a relatively large a, namely a = 10. Figure 1.4 shows the result of 
this calculation, which appears quite reasonable. 




Figure 1.4: Solution of (1.5) with initial data (1.9) using (1.7) with 100 points in 
the space grid and a = 10. To avoid an over-dense graph not all the points in the 
numerical grid are plotted (enough points to show all the relevant details are kept). 



In the second calculation, we take a smaller value a = 6. This makes the corners more substantial 
(though still pretty weak). Figure 1.5 shows the result of this last calculation, which is now not 
reasonable at all. It is quite clear that, just as in the prior example, the small errors that are 
triggered by the corners are amplified by the scheme (so we observe grid scale oscillations near 
x = ±1 towards the end of the run). 

Finally, we point out that, if the calculations are run for times longer than < t < 0.5, even the one 
with a = 10 eventually shows grid scale oscillations. These grow exponentially in time and pretty 
soon dominate the whole solution (not just the neighborhood of x = ±1) with huge amplitudes. 
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Figure 1.5: Solution of (1.5) with initial data (1.9) using (1.7) with 100 points in 
the space grid and a = 6. To avoid an over-dense graph not all the points in the 
numerical grid are plotted (enough points to show all the relevant details are kept). 



The next section gives a detailed explanation of why this is happening. 

2 von Neumann stability analysis for PDE's. 

In this section we introduce the von Neumann stability analysis technique, that can be used to 
analyze numerical schemes and predict when the behavior observed in the prior section will occur. 
There are two basic concepts useful in understanding numerical schemes. These are the notions of 
consistency and stability. For a numerical scheme to be useful it must be both consistent and 
stable. It is very important to realize that these two notions are independent. 

Consistency simply means that, as Ax and At vanish, the solutions of the equation must satisfy 
the numerical scheme with errors that vanish. This is in fact what equation (1.6) tells us about 
the scheme in (1.7). Consistency guarantees that the scheme truly approximates the equation we 
intend to solve with it (and not something else). 
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Stability simply means that the scheme does not amplify errors. Obviously this is very important, 
since errors are impossible to avoid in any numerical calculation. In fact, even in the ideal case 
of infinite precision, we still have to deal with discretization errors — i.e. the O terms in (1.6). 
Clearly, if errors are amplified, pretty soon they will dominate any computation (making it useless). 

As it turns out, for linear constant coefficient schemes such as (1.7), a complete stability 
analysis is possible, because the numerical algorithm equations can be solved exactly by separa- 
tion of variables. This means then that any solution of the scheme can be written as a superposition 
of Fourier modes. These Fourier modes are solutions of the form 

u{ = U G j e ikn and v J n = V G j e ikn , (2.1) 

where U, V, G and k are constants (with k real). Generally double sequences like this will be solu- 
tions provided G, U and V are restricted by some functional relations of the form G = G(k, Ax, At), 
U = U(k, Ax, At) and V = V(k, Ax, At) — below we carry through the calculations for the specific 
example of (1.7). 

G is called the Growth Factor. It is clear that: 

' (2-2) 



for stability ||C|| < 1 is needed for all k. 



Else some modes will be amplified by a factor G in each time step, eventually dominating the 
solution. A scheme is called stable if the stability condition ||G|| < 1 can be satisfied with (perhaps) 
a restriction on the time step of the form < At < t(Ax), where r is a positive function of its 
argument. Notice that restrictions of this latter form allow arbitrarily small time and space steps, 
which are needed to be able to compute the solution with any required degree of accuracy (how 
small is determined by how well consistency is satisfied, which determines the size of the errors for 
any given At and Ax). 

Remark 2.1 The parameter k is the wavenumber of the mode, related to the wavelength A in 
space 1 by A = {2'wAx)/k. For the particular case of periodic problems (such as the ones consid- 
ered in examples 1.1 and 1.2), the Fourier modes (2.1) must also satisfy the periodicity condition. 
That is, one must have A = T/t, where I is an integer and T is the period in space. Since in this 
case one would normally take Ax = T/N , where N is a large natural number, the acceptable values 
for k end up restricted to the set 

k = k t = * X 1= -^l and \ = \ t = -, with 0<£<N. (2.3) 

Here the upper bound N on I follows from the fact that ki and ki +N give the same Fourier mode in 
(2.1); thus there is no reason to keep both. 

We note that (due to the fact that the numerical scheme only samples the solution at a discrete set 
{ x n} of points in space) there is a certain trickiness in the interpretation of the wavelengths 
Xi above. Clearly, I = corresponds to a solution independent of x and l= \ corresponds to the 
fundamental mode with wavelength T in x. As I continues to increase harmonics of this fundamental 
mode appear, with wavelengths T/2, T/3 . . . However, this process cannot continue forever, since 



k 

1 Write the argument kn in the exponentials in (2.1) as kn = — — (x n — xq), using (1.3). 

Ax 
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the numerical grid cannot resolve arbitrarily small wavelengths. In fact, the shortest wavelength that 
can be resolved corresponds to I = N/2 with Xi = 2 Ax (grid size oscillations, with period 2 in n: the 
solution alternates between two values on the grid). To see this recall that ki and ki + N give the same 
Fourier mode in (2.1). Thus the mode (N — I) has the same wavelength as the mode —I, i.e. T jl. 
This means that, after I = N/2 the wavelengths start increasing, to reach back the fundamental 
mode at I = N — 1. Each wavelength then actually appears twice in the range 1 < I < N. 

We should not be too surprised by the fact that each wavelength appears twice in the range 1 < I < N. 
Notice that the modes in (2.1) are complex valued (except when k is a multiple of 2tt). Thus, to be 
real valued any solution should include both the modes and their complex conjugates. However, the 
mode conjugate to the one with k = kt above in (2.3) is the mode with k = k-t, which is precisely 
the same as the mode with k = k^-i- 

In any numerical calculation it is the modes with wavelengths of the order of the grid size Ax 
(i.e. I close to N/2) that are worrisome in terms of instabilities. These modes cannot be expected 
to represent accurately any true feature of the real solution one is trying to compute 2 and should 
not have any significant presence in the numerical solution. Thus, it is very important that they 
not be amplified by the scheme. In fact, generally it is desirable to have them damped, since they 
mostly represent numerical "noise" generated by all the approximations implicit in any numerical 
calculation. 

On the other hand, the modes with wavelengths much bigger than Ax (that is, I ~ or I ~ N in 
(2.3)) should be treated "accurately" by the scheme. By this we mean that their time evolution 
(given by the factors G 3 in (2.1)) should be as close as possible to the one provided by the PDE the 
scheme approximates. This is what consistency is all about. 

Consider now the special case of the algorithm (1.7). To see under which conditions (2.1) 
is a solution, substitute this form into (1.7). Dividing by the common factor G 3 e lkn it follows that 



Clearly an eigenvalue equation AY = GY, with eigenvalue G, eigenvector Y = (U, V) T and matrix 
of coefficients 



GU = U + AtV and GV = V + 



At 



(e ik -2 + e~ ik )U. 



(Ax) 2 




From the characteristic equation det(A — G) = 0, then 




(2.4) 




(2.5) 



which is always bigger than one. 



2 Recall (1.4), which makes sense in terms of approximating the solution only if Ax is much smaller than any 
distance over which the solution changes significantly. 
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Notice that the maximum amplification for the scheme (1.7) occurs — as follows from (2.5) 
— for k = 7r. This corresponds to I = N/2 in (2.3), i.e.: grid size oscillations with A = 2 Ax. 
In this case 

\\G\\=G M = yJl + 4v, (2-6) 



where 



n = (At/ Ax) 2 . 



For (1.7), the amplitude of the grid size oscillations grows like G 3 M . Thus 
we can write for the amplification factor A 2 = A 2 (t) (for the period 2 Ax mode) 

A 2 = exp(i^) , (2.7) 

where we have used j = t/ At. In particular (in examples 1.1 and 1.2 earlier) we took Ax = 2 At 
and At = 1/N, so that 

. In 2 . jv t , 
A 2 = exp(— Nt) = 2—. (2.8) 

We will now use these results to explain the behavior observed earlier in figures 1.1 through 1.5. 

Remark 2.2 Consider first example 1.1, with the initial data for scheme (1.7) given by 

1 ( , 277.7T \ n 

Un = 2 V C ° S (~/\P ) Vn = ' 

These data correspond to a superposition of just three modes in (2.1), with k = k , k = k\ and 
k = ~ k N _i in (2.3). Thus, the exact solution for the scheme equations is rather simple 
and has the form 

9 ] + 9 j . 2n7^ ^^ . j 9 j ~ 9 j . ,2nyr . tt 

— 2 — C0S (^F') n= 2i v C0S (^F' ' /or ^ = 1 + ?sin (^)' ( 2 - 9 ) 

where v is a constant and g denotes the complex conjugate of g. Of course, g and g are the values 
G in (2.4) takes for k = h = 2ir/N. 

Notice that the exact solution (2.9) does not exhibit any catastrophic growth of grid size oscillations, 
as was observed in example 1.1. However, the results displayed in figures 1.1 through 1.3 do not 
correspond to the exact solution above but to actual computations using the scheme in (1.7) — which 
were done using double precision floating point arithmetic (MatLab's default). The round off errors 
introduced by the finite precision of the calculations introduces (very small) perturbations into the 
exact solution above, which the scheme then evolves in time just as if they were part of the solution. 

To understand what the scheme does with the perturbations introduced by the finite precision, de- 
compose them into a sum over the modes in (2.1). This sum will generally include all the modes, 
in particular the highly amplified ones with grid size wavelengths. Consider then what would happen 
with the solution of the scheme if we add to the initial data above 3 a small amount of the component 
corresponding to the maximum amplification rate above in (2.6). Let the amplitude of this compo- 
nent be e, where e has (roughly) the size of the expected errors. Actually, e should be a little smaller 
than the round off errors that occur, since not all the errors get projected into the fastest growing 

as a good ballpark figure for the calculations in section 1 




modes. Thus take 



e = O(10~ 17 ) 



and use (2.8) above to explain the behavior observed in figures 1.1 through 1.3, as follows: 



3 Which has only components corresponding to I = 0, I = 1 and I = N — 1 in (2.3). 
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1. First, for N = 40, (2.8) gives A 2 ~ 1.1 x 10 12 for the final time t = 2. This is not enough to 
compensate for the smallness of e and the numerical solution is well described by (2.9). 

Notice that (2.9) is not periodic in time; since the wave amplitude in u behaves like Re(gi), 
which grows as j grows. In fact, 2 N = 80 steps are needed to reach the final time t = 2 and 
it is easy to check that 

Re(g m ) = Re j (l + i M^)) ^ J ~ 1-28 . 

This agrees quite well with the ~ 30% growth in the wave amplitude observed in figure 1.1. 

2. Second, for N = 57, (2.8) gives A 2 ~ 1.4 x 10 17 for the final time t = 2. This is about the 
same as e -1 and agrees with the fact that grid oscillations of 0(1) amplitude are observed in 
figure 1.2. 

3. Third, for N = 80, (2.8) gives A 2 « 1.2 x 10 24 for the final time t = 2. This is about 10 7 
times bigger than e -1 , which (again) agrees pretty well with the observed amplitude of the grid 
size oscillations in figure 1.3. 

4- Finally, it is not just the mode with I = N/2 in (2.3) that gets a large amplification factor by 
the scheme. All the ones with I ~ N/2 do and should thus be present in the solution. It is 
well known that when sinusoidals with close wavenumbers are added, "beats" with wavenumbers 
equal to the difference in wavenumbers occur. Thus, in this case we should observe "beats" with 
wavenumbers low multiples of k\ = 2tt/N — which, indeed, are quite obvious in figure 1.3. 

Remark 2.3 Now consider example 1.2, w here N = 100 and < t < 0.5. Then, for the time 
t = 0.5, equation (2.8) gives A 2 ~ 3.4 x 10' 



In this case the initial data has components in all the modes < I < N in (2.3). In fact, be- 
cause of the corners at x = ±1, the amplitude present in the higher modes is relatively large. The 
strength of these corners can be measured by the jump in the derivative of the initial data there: 



J(a) = 4a ln(10) 10 a . For moderate* size a, J(a) pretty much determines how much amplitude 

there is in the higher modes. Now J(10) ~ 9.2 x 10~ 9 and J(6) ~ 5.5 x 10~ 5 . Thus, from the value 
of A 2 above, it should be clear why in figure 1.4 (corresponding to a = 10) the solution exhibits no 
detectable oscillations, while in figure 1.5 (corresponding to a = 6) they show up. 

Notice that in this case it is also true that it is not just the mode with I = N/2 in (2.3) that gets a 
large amplification factor by the scheme. All the neighboring ones are also present. However, now 
their amplitudes and phases are all correlated because they (mostly) are generated by the corner in 
the initial data. Thus they interfere with each other in ways subtler than the mere beating observed in 
the prior example; i.e.: the pattern of grid size oscillations has a clear maximum near the positions 
of the corners in figure 1.5. 

In the next section we will discuss a simple strategy to stabilize numerical schemes, to get rid 
of numerical oscillations and other undesirable effects. The strategy is based on the introduction 
of artificial (numerical) dissipation to (selectively) damp the higher modes, without significantly 
affecting the lower modes (where a consistent scheme should behave properly — see remark 2.4). 



4 When a is large, the corner is very weak and the dominant contribution to the mode amplitudes comes from the 
smooth part of the initial data (which yields very little amplitude in the high modes). 
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Remark 2.4 Finally, going back now to the last paragraph in remark 2.1, consider the behavior of 
G in (2.4) for k small. Namely 

G = 1±< £* + ° (£;*•)• (2 - 10) 

This should be compared with the behavior of the exact solution for the wave equation (1.1) — see 
remark 1.1 — which evolves Fourier modes according to the rule 



u oc 



exp J i^( x n ± tj) J oc exp ji (kn ± ^kj^j J . 



Thus the exact evolution corresponds to a factor G given by 

G exact = exp (±^) = l±i^k + ((|£*) 2 ) • (2.11) 

This should be compared with (2.10) above. It is clear then that (for k small) G is correct up to 
small terms in k, which is an alternative way of verifying that the scheme (1.7) is consistent. 

3 Numerical Viscosity and Stabilized Scheme. 

FILL IN HERE THE GOOD SCHEME EQUATIONS. (3.1) 

Notation used for Good Scheme in MatLab: r\ = (At/ Ax) 2 and v = At/ Ax 2 . 
Next the figures that go with the good scheme. 

4 Reference. 

For more information regarding stability of numerical schemes (and many other useful numerical 
topics) a good all-around practical reference is Numerical Recipes, The Art of Scientific Computing 
by W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery. Cambridge U. Press, New 
York, 1992. 
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Figure 3.1: Solution of (1.5) with initial data (1.8) using the corrected scheme (3.1) 
with 55 points in the space grid. To avoid an over-dense graph not all the points 
in the numerical grid are plotted. However, enough points to show all the relevant 
details are kept. 
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Figure 3.2: Solution of (1.5) with initial data (1.8) using the corrected scheme (3.1) 
with 190 points in the space grid. To avoid an over-dense graph not all the points 
in the numerical grid are plotted. However, enough points to show all the relevant 
details are kept. 



